
library(mlogit)
library(gmnl)
library(aspace)


###### table 1 ######
df <- read.csv("df_wide.csv")
df <-df[is.na(df$S1)==F,]
colnames(df)
df$gender <- df$Q40
df$race <- df$Q41
df$marriage <- df$Q45
df$urban <- df$Q44
df$adult <- df$Q46
df$child <- df$Q47
df$political <- df$Q87


unique(df$Q42[!grepl("^\\s*\\d+(\\.\\d+)?\\s*$", df$Q42)])
df[df$Q42 == ",74",]$Q42  <- "74"
df$age <- as.numeric(df$Q42)

df$edu <- df$Q48
df$income <- df$Q49


cat_vars <- c("gender", "race", "edu", "income", "marriage", "urban", "adult", "child", "political")
cont_vars <- c("age")

var <- c("gender", "race", "edu", "income", "age")

summary_cat(vars = cat_vars,data = df)
chifun(vars = cat_vars, treatment = "Treat", data = df)


summary_cont(vars = cont_vars, data = df)
aovfun(vars = cont_vars, treatment = "Treat", data = df)

summary_cont(vars = var, data = df)




df$S1_self <- 85
df$S2_self <- ifelse(df$S2 == 1, 85, ifelse(df$S2 == 2, 87, 
                                            ifelse(df$S2 == 3, 89, ifelse(df$S2 == 4, 91,
                                                                          ifelse(df$S2 == 5, 93, ifelse(df$S2 == 6, 94,
                                                                                                        ifelse(df$S2 == 7, 96, ifelse(df$S2 == 8, 98,
                                                                                                                                      ifelse(df$S2 == 9, 100, -999)))))))))
table(df$S2_self)
df$S3_self <- ifelse(df$S3 == 1, 50, ifelse(df$S3 == 2, 54, 
                                            ifelse(df$S3 == 3, 59, ifelse(df$S3 == 4, 63,
                                                                          ifelse(df$S3 == 5, 68, ifelse(df$S3 == 6, 72,
                                                                                                        ifelse(df$S3 == 7, 76, ifelse(df$S3 == 8, 81,
                                                                                                                                      ifelse(df$S3 == 9, 85, -999)))))))))
df$S4_self <- ifelse(df$S4 == 1, 50, ifelse(df$S4 == 2, 54, 
                                            ifelse(df$S4 == 3, 59, ifelse(df$S4 == 4, 63,
                                                                          ifelse(df$S4 == 5, 68, ifelse(df$S4 == 6, 72,
                                                                                                        ifelse(df$S4 == 7, 76, ifelse(df$S4 == 8, 81,
                                                                                                                                      ifelse(df$S4 == 9, 85, -999)))))))))
df$S5_self <- ifelse(df$S5 == 1, 100, ifelse(df$S5 == 2, 94, 
                                             ifelse(df$S5 == 3, 88, ifelse(df$S5 == 4, 81,
                                                                           ifelse(df$S5 == 5, 75, ifelse(df$S5 == 6, 69,
                                                                                                         ifelse(df$S5 == 7, 63, ifelse(df$S5 == 8, 56,
                                                                                                                                       ifelse(df$S5 == 9, 50, -999)))))))))
df$S6_self <- ifelse(df$S6 == 1, 100, ifelse(df$S6 == 2, 98, 
                                             ifelse(df$S6 == 3, 96, ifelse(df$S6 == 4, 94,
                                                                           ifelse(df$S6 == 5, 93, ifelse(df$S6 == 6, 91,
                                                                                                         ifelse(df$S6 == 7, 89, ifelse(df$S6 == 8, 87,
                                                                                                                                       ifelse(df$S6 == 9, 85, -999)))))))))
df$S1_other <- ifelse(df$S1 == 1, 85, ifelse(df$S1 == 2, 76, 
                                             ifelse(df$S1 == 3, 68, ifelse(df$S1 == 4, 59,
                                                                           ifelse(df$S1 == 5, 50, ifelse(df$S1 == 6, 41,
                                                                                                         ifelse(df$S1 == 7, 33, ifelse(df$S1 == 8, 24,
                                                                                                                                       ifelse(df$S1 == 9, 15, -999)))))))))
df$S2_other <- ifelse(df$S2 == 1, 15, ifelse(df$S2 == 2, 19, 
                                             ifelse(df$S2 == 3, 24, ifelse(df$S2 == 4, 28,
                                                                           ifelse(df$S2 == 5, 33, ifelse(df$S2 == 6, 37,
                                                                                                         ifelse(df$S2 == 7, 41, ifelse(df$S2 == 8, 46,
                                                                                                                                       ifelse(df$S2 == 9, 50, -999)))))))))
df$S3_other <- ifelse(df$S3 == 1, 100, ifelse(df$S3 == 2, 98, 
                                              ifelse(df$S3 == 3, 96, ifelse(df$S3 == 4, 94,
                                                                            ifelse(df$S3 == 5, 93, ifelse(df$S3 == 6, 91,
                                                                                                          ifelse(df$S3 == 7, 89, ifelse(df$S3 == 8, 87,
                                                                                                                                        ifelse(df$S3 == 9, 85, -999)))))))))
df$S4_other <- ifelse(df$S4 == 1, 100, ifelse(df$S4 == 2, 89, 
                                              ifelse(df$S4 == 3, 79, ifelse(df$S4 == 4, 68,
                                                                            ifelse(df$S4 == 5, 58, ifelse(df$S4 == 6, 47,
                                                                                                          ifelse(df$S4 == 7, 36, ifelse(df$S4 == 8, 26,
                                                                                                                                        ifelse(df$S4 == 9, 15, -999)))))))))
df$S5_other <- ifelse(df$S5 == 1, 50, ifelse(df$S5 == 2, 56, 
                                             ifelse(df$S5 == 3, 63, ifelse(df$S5 == 4, 69,
                                                                           ifelse(df$S5 == 5, 75, ifelse(df$S5 == 6, 81,
                                                                                                         ifelse(df$S5 == 7, 88, ifelse(df$S5 == 8, 94,
                                                                                                                                       ifelse(df$S5 == 9, 100, -999)))))))))
df$S6_other <- ifelse(df$S6 == 1, 50, ifelse(df$S6 == 2, 54, 
                                             ifelse(df$S6 == 3, 59, ifelse(df$S6 == 4, 63,
                                                                           ifelse(df$S6 == 5, 68, ifelse(df$S6 == 6, 72,
                                                                                                         ifelse(df$S6 == 7, 76, ifelse(df$S6 == 8, 81,
                                                                                                                                       ifelse(df$S6 == 9, 85, -999)))))))))
df$self_average <- (df$S1_self + df$S2_self + df$S3_self + df$S4_self + df$S5_self + df$S6_self)/6

df$other_average <- (df$S1_other + df$S2_other + df$S3_other + df$S4_other + df$S5_other + df$S6_other)/6


df$svo <- atan_d((df$other_average - 50)/(df$self_average - 50))

df$social_cate <- ifelse(df$svo >57.15, "Altruism", ifelse(df$svo >22.45, "Prosociality",
                                                           ifelse(df$svo > -12.04, "Individualism",
                                                                  ifelse(df$svo <= -12.04, "Competitiveness", "Not_find"))))

df$indi_comp <- ifelse(df$social_cate == "Individualism" | df$social_cate == "Competitiveness", 1,0)

df$indi <- ifelse(df$social_cate == "Individualism", 1,0)

df$comp <- ifelse(df$social_cate == "Competitiveness", 1,0)

df$waste <- df$Q27_1
df$indi_comp <- as.factor(df$indi_comp)
# Compute the analysis of variance
res.aov <- aov(waste ~ indi_comp, data = df)
# Summary of the analysis
summary(res.aov)

TukeyHSD(res.aov)

summary(df[df$indi_comp==1,]$waste)

summary(df$waste)

####### analysis #######
df.long <- read.csv("df_long.csv")

df.long <-df.long[is.na(df.long$S1)==F,]

df.long$nobuy <- ifelse(df.long$alt == 3, 1,0)


df.long$NEP_new_1 <- ifelse(df.long$Q29_1 == 1, 5, 
                              ifelse(df.long$Q29_1 == 2, 4,
                                     ifelse(df.long$Q29_1 == 3, 3, 
                                            ifelse(df.long$Q29_1 == 4, 2,1))))
df.long$NEP_new_2 <- df.long$Q29_2
df.long$NEP_new_3 <- ifelse(df.long$Q29_3 == 1, 5, 
                              ifelse(df.long$Q29_3 == 2, 4,
                                     ifelse(df.long$Q29_3 == 3, 3, 
                                            ifelse(df.long$Q29_3 == 4, 2,1))))
df.long$NEP_new_4 <- df.long$Q29_4
df.long$NEP_new_5 <- ifelse(df.long$Q29_5 == 1, 5, 
                              ifelse(df.long$Q29_5 == 2, 4,
                                     ifelse(df.long$Q29_5 == 3, 3, 
                                            ifelse(df.long$Q29_5 == 4, 2,1))))


df.long$NEP_new_6 <- df.long$Q30_1
df.long$NEP_new_7 <- ifelse(df.long$Q30_2 == 1, 5, 
                              ifelse(df.long$Q30_2 == 2, 4,
                                     ifelse(df.long$Q30_2 == 3, 3, 
                                            ifelse(df.long$Q30_2 == 4, 2,1))))
df.long$NEP_new_8 <- df.long$Q30_3
df.long$NEP_new_9 <- ifelse(df.long$Q30_4 == 1, 5, 
                              ifelse(df.long$Q30_4 == 2, 4,
                                     ifelse(df.long$Q30_4 == 3, 3, 
                                            ifelse(df.long$Q30_4 == 4, 2,1))))
df.long$NEP_new_10 <- df.long$Q30_5



df.long$NEP_new_11 <- ifelse(df.long$Q31_1 == 1, 5, 
                               ifelse(df.long$Q31_1 == 2, 4,
                                      ifelse(df.long$Q31_1 == 3, 3, 
                                             ifelse(df.long$Q31_1 == 4, 2,1))))

df.long$NEP_new_12 <- df.long$Q31_2
df.long$NEP_new_13 <- ifelse(df.long$Q31_3 == 1, 5, 
                               ifelse(df.long$Q31_3 == 2, 4,
                                      ifelse(df.long$Q31_3 == 3, 3, 
                                             ifelse(df.long$Q31_3 == 4, 2,1))))
df.long$NEP_new_14 <- df.long$Q31_4
df.long$NEP_new_15 <- ifelse(df.long$Q31_5 == 1, 5, 
                               ifelse(df.long$Q31_5 == 2, 4,
                                      ifelse(df.long$Q31_5 == 3, 3, 
                                             ifelse(df.long$Q31_5 == 4, 2,1))))

df.long$NEP_total <- df.long$NEP_new_1 + df.long$NEP_new_2 + df.long$NEP_new_3 + df.long$NEP_new_4 +
  df.long$NEP_new_5 + df.long$NEP_new_6 + df.long$NEP_new_7 + df.long$NEP_new_8 +
  df.long$NEP_new_9 + df.long$NEP_new_10 + df.long$NEP_new_11 + df.long$NEP_new_12 + 
  df.long$NEP_new_13 + df.long$NEP_new_14 + df.long$NEP_new_15

median(df.long$NEP_total)
mean(df.long$NEP_total)

nrow(df.long[df.long$NEP_total < median(df.long$NEP_total),])
nrow(df.long[df.long$NEP_total >= median(df.long$NEP_total),])

nrow(df.long[df.long$NEP_total < mean(df.long$NEP_total),])
nrow(df.long[df.long$NEP_total >= mean(df.long$NEP_total),])

# use >= median as high NEP

df.long$hig_NEP <- ifelse(df.long$NEP_total >= median(df.long$NEP_total), 1,0)

sum(df.long$hig_NEP)


df.long$S1_self <- 85
df.long$S2_self <- ifelse(df.long$S2 == 1, 85, ifelse(df.long$S2 == 2, 87, 
                                                      ifelse(df.long$S2 == 3, 89, ifelse(df.long$S2 == 4, 91,
                                                                                         ifelse(df.long$S2 == 5, 93, ifelse(df.long$S2 == 6, 94,
                                                                                                                            ifelse(df.long$S2 == 7, 96, ifelse(df.long$S2 == 8, 98,
                                                                                                                                                               ifelse(df.long$S2 == 9, 100, -999)))))))))
table(df.long$S2_self)
df.long$S3_self <- ifelse(df.long$S3 == 1, 50, ifelse(df.long$S3 == 2, 54, 
                                                      ifelse(df.long$S3 == 3, 59, ifelse(df.long$S3 == 4, 63,
                                                                                         ifelse(df.long$S3 == 5, 68, ifelse(df.long$S3 == 6, 72,
                                                                                                                            ifelse(df.long$S3 == 7, 76, ifelse(df.long$S3 == 8, 81,
                                                                                                                                                               ifelse(df.long$S3 == 9, 85, -999)))))))))
df.long$S4_self <- ifelse(df.long$S4 == 1, 50, ifelse(df.long$S4 == 2, 54, 
                                                      ifelse(df.long$S4 == 3, 59, ifelse(df.long$S4 == 4, 63,
                                                                                         ifelse(df.long$S4 == 5, 68, ifelse(df.long$S4 == 6, 72,
                                                                                                                            ifelse(df.long$S4 == 7, 76, ifelse(df.long$S4 == 8, 81,
                                                                                                                                                               ifelse(df.long$S4 == 9, 85, -999)))))))))
df.long$S5_self <- ifelse(df.long$S5 == 1, 100, ifelse(df.long$S5 == 2, 94, 
                                                      ifelse(df.long$S5 == 3, 88, ifelse(df.long$S5 == 4, 81,
                                                                                         ifelse(df.long$S5 == 5, 75, ifelse(df.long$S5 == 6, 69,
                                                                                                                            ifelse(df.long$S5 == 7, 63, ifelse(df.long$S5 == 8, 56,
                                                                                                                                                               ifelse(df.long$S5 == 9, 50, -999)))))))))
df.long$S6_self <- ifelse(df.long$S6 == 1, 100, ifelse(df.long$S6 == 2, 98, 
                                                      ifelse(df.long$S6 == 3, 96, ifelse(df.long$S6 == 4, 94,
                                                                                         ifelse(df.long$S6 == 5, 93, ifelse(df.long$S6 == 6, 91,
                                                                                                                            ifelse(df.long$S6 == 7, 89, ifelse(df.long$S6 == 8, 87,
                                                                                                                                                               ifelse(df.long$S6 == 9, 85, -999)))))))))
df.long$S1_other <- ifelse(df.long$S1 == 1, 85, ifelse(df.long$S1 == 2, 76, 
                                                      ifelse(df.long$S1 == 3, 68, ifelse(df.long$S1 == 4, 59,
                                                                                         ifelse(df.long$S1 == 5, 50, ifelse(df.long$S1 == 6, 41,
                                                                                                                            ifelse(df.long$S1 == 7, 33, ifelse(df.long$S1 == 8, 24,
                                                                                                                                                               ifelse(df.long$S1 == 9, 15, -999)))))))))
df.long$S2_other <- ifelse(df.long$S2 == 1, 15, ifelse(df.long$S2 == 2, 19, 
                                                      ifelse(df.long$S2 == 3, 24, ifelse(df.long$S2 == 4, 28,
                                                                                         ifelse(df.long$S2 == 5, 33, ifelse(df.long$S2 == 6, 37,
                                                                                                                            ifelse(df.long$S2 == 7, 41, ifelse(df.long$S2 == 8, 46,
                                                                                                                                                               ifelse(df.long$S2 == 9, 50, -999)))))))))
df.long$S3_other <- ifelse(df.long$S3 == 1, 100, ifelse(df.long$S3 == 2, 98, 
                                                      ifelse(df.long$S3 == 3, 96, ifelse(df.long$S3 == 4, 94,
                                                                                         ifelse(df.long$S3 == 5, 93, ifelse(df.long$S3 == 6, 91,
                                                                                                                            ifelse(df.long$S3 == 7, 89, ifelse(df.long$S3 == 8, 87,
                                                                                                                                                               ifelse(df.long$S3 == 9, 85, -999)))))))))
df.long$S4_other <- ifelse(df.long$S4 == 1, 100, ifelse(df.long$S4 == 2, 89, 
                                                      ifelse(df.long$S4 == 3, 79, ifelse(df.long$S4 == 4, 68,
                                                                                         ifelse(df.long$S4 == 5, 58, ifelse(df.long$S4 == 6, 47,
                                                                                                                            ifelse(df.long$S4 == 7, 36, ifelse(df.long$S4 == 8, 26,
                                                                                                                                                               ifelse(df.long$S4 == 9, 15, -999)))))))))
df.long$S5_other <- ifelse(df.long$S5 == 1, 50, ifelse(df.long$S5 == 2, 56, 
                                                       ifelse(df.long$S5 == 3, 63, ifelse(df.long$S5 == 4, 69,
                                                                                          ifelse(df.long$S5 == 5, 75, ifelse(df.long$S5 == 6, 81,
                                                                                                                             ifelse(df.long$S5 == 7, 88, ifelse(df.long$S5 == 8, 94,
                                                                                                                                                                ifelse(df.long$S5 == 9, 100, -999)))))))))
df.long$S6_other <- ifelse(df.long$S6 == 1, 50, ifelse(df.long$S6 == 2, 54, 
                                                       ifelse(df.long$S6 == 3, 59, ifelse(df.long$S6 == 4, 63,
                                                                                          ifelse(df.long$S6 == 5, 68, ifelse(df.long$S6 == 6, 72,
                                                                                                                             ifelse(df.long$S6 == 7, 76, ifelse(df.long$S6 == 8, 81,
                                                                                                                                                                ifelse(df.long$S6 == 9, 85, -999)))))))))
df.long$self_average <- (df.long$S1_self + df.long$S2_self + df.long$S3_self + df.long$S4_self + df.long$S5_self + df.long$S6_self)/6

df.long$other_average <- (df.long$S1_other + df.long$S2_other + df.long$S3_other + df.long$S4_other + df.long$S5_other + df.long$S6_other)/6


df.long$svo <- atan_d((df.long$other_average - 50)/(df.long$self_average - 50))

df.long$social_cate <- ifelse(df.long$svo >57.15, "Altruism", ifelse(df.long$svo >22.45, "Prosociality",
                                                                     ifelse(df.long$svo > -12.04, "Individualism",
                                                                            ifelse(df.long$svo <= -12.04, "Competitiveness", "Not_find"))))

df.long$indi_comp <- ifelse(df.long$social_cate == "Individualism" | df.long$social_cate == "Competitiveness", 1,0)

df.long$indi <- ifelse(df.long$social_cate == "Individualism", 1,0)

df.long$comp <- ifelse(df.long$social_cate == "Competitiveness", 1,0)

summary(df.long$Q27_1)

nrow(df.long[df.long$Q27_1 < median(df.long$Q27_1),])

df.long$low_waste <- ifelse(df.long$Q27_1 < median(df.long$Q27_1), 1,0)

100*table(df.long[df.long$low_waste==1,]$social_cate)/sum(table(df.long[df.long$low_waste==1,]$social_cate))
100*table(df.long[df.long$low_waste==0,]$social_cate)/sum(table(df.long[df.long$low_waste==0,]$social_cate))


100*table(df.long[df.long$Q27_1 < mean(df.long$Q27_1),]$social_cate)/sum(table(df.long[df.long$Q27_1 < mean(df.long$Q27_1),]$social_cate))
100*table(df.long[df.long$Q27_1 > mean(df.long$Q27_1),]$social_cate)/sum(table(df.long[df.long$Q27_1 > mean(df.long$Q27_1),]$social_cate))


table(df.long$social_cate)

n.ml<-mlogit.data(data = df.long, 
                  choice ="Choice", 
                  shape = "long", 
                  alt.var = "alt", 
                  reflevel = 3,
                  id.var = "ID"
)

ml.wtp <-mlogit.data(data = df.long, 
                  choice ="Choice", 
                  shape = "long", 
                  alt.var = "alt", 
                  reflevel = 3,
                  id.var = "ID",
                  opposite = "price"
)

table(df.long$Treat)

####### mnl  ########
mnl.control <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
                    reflevel = 3,
                    model = "mnl",  
                    print.init = TRUE,
                    print.level = 2,
                    data = n.ml,
                    subset = n.ml$Treat == 0
                    
)
summary(mnl.control)

mnl.t1 <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
               reflevel = 3,
               model = "mnl",  
               print.init = TRUE,
               print.level = 2,
               data = n.ml,
               subset = n.ml$Treat == 1
               
)
summary(mnl.t1)

mnl.t2 <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
               reflevel = 3,
               model = "mnl",  
               print.init = TRUE,
               print.level = 2,
               data = n.ml,
               subset = n.ml$Treat == 2
               
)
summary(mnl.t2)



####### Table 2 #######
wtp.control <-coef(mnl.control)/(0-coef(mnl.control)[2])

start.control <- c(1, wtp.control[c(1,3,4)],0,rep(0.1,3),0.1,0)
wtp.rpl.control <- gmnl(formula = Choice ~ price + nobuy + gene  + color|0|0|0|1,
                              reflevel = 3,
                              model="gmnl",
                              data=ml.wtp,
                              subset = ml.wtp$Treat==0,
                              R=500,
                              ranp = c(gene = "n", color = "n"),
                              correlation = T,
                              fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                              panel = T, 
                              haltons = NA,
                              start = start.control,
                              print.init = 2,
                              print.level =2,
                              method = "bhhh",
                              rterlim = 500
                              
)
summary(wtp.rpl.control)
AIC(wtp.rpl.control)
BIC(wtp.rpl.control)
save(wtp.rpl.control, file = "wtp.rpl.control.Rdata")
load("wtp.rpl.control.Rdata")

se.cov.gmnl(wtp.rpl.control,sd=T)


wtp.t1 <-coef(mnl.t1)/(0-coef(mnl.t1)[2])

start.t1 <- c(1, wtp.t1[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t1 <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                         reflevel = 3,
                         model="gmnl",
                         data=ml.wtp,
                         subset = ml.wtp$Treat==1, 
                         R=500,
                         ranp = c(gene = "n", color = "n"),
                         correlation = T,
                         fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                         panel = T, 
                         haltons = NA,
                         start = start.t1,
                         print.init = 2,
                         print.level =2,
                         method = "bhhh",
                         rterlim = 500
                         
)
summary(wtp.rpl.t1)
AIC(wtp.rpl.t1)
BIC(wtp.rpl.t1)
save(wtp.rpl.t1, file = "wtp.rpl.t1.Rdata")
load("wtp.rpl.t1.Rdata")

se.cov.gmnl(wtp.rpl.t1,sd=T)

wtp.t2 <-coef(mnl.t2)/(0-coef(mnl.t2)[2])

start.t2 <- c(1, wtp.t2[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t2 <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                         reflevel = 3,
                         model="gmnl",
                         data=ml.wtp,
                         subset = ml.wtp$Treat==2,
                         R=500,
                         ranp = c(gene = "n", color = "n"),
                         correlation = T,
                         fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                         panel = T, 
                         haltons = NA,
                         start = start.t2,
                         print.init = 2,
                         print.level =2,
                         method = "bhhh",
                         rterlim = 500
                         
)
summary(wtp.rpl.t2)
AIC(wtp.rpl.t2)
BIC(wtp.rpl.t2)
save(wtp.rpl.t2, file = "wtp.rpl.t2.Rdata")
load("wtp.rpl.t2.Rdata")

se.cov.gmnl(wtp.rpl.t2,sd=T)

######## table 3 #######


model1 <- c("wtp.rpl.control","wtp.rpl.control","wtp.rpl.t1")
model2 <- c("wtp.rpl.t1","wtp.rpl.t2","wtp.rpl.t2")
vari_compair <- c("gene", "color")
kr.poe(model1, model2, vari_compair)


####### subsample analysis #######
#### indi comp ####
mnl.control.indi_comp <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
                    reflevel = 3,
                    model = "mnl",  
                    print.init = TRUE,
                    print.level = 2,
                    data = n.ml,
                    subset = n.ml$Treat == 0 & n.ml$indi_comp ==1
                    
)
summary(mnl.control.indi_comp)

mnl.control.other <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
                              reflevel = 3,
                              model = "mnl",  
                              print.init = TRUE,
                              print.level = 2,
                              data = n.ml,
                              subset = n.ml$Treat == 0 & n.ml$indi_comp ==0
                              
)
summary(mnl.control.other)


mnl.t1.indi_comp <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
               reflevel = 3,
               model = "mnl",  
               print.init = TRUE,
               print.level = 2,
               data = n.ml,
               subset = n.ml$Treat == 1 & n.ml$indi_comp ==1
               
)
summary(mnl.t1.indi_comp)

mnl.t1.other <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
                         reflevel = 3,
                         model = "mnl",  
                         print.init = TRUE,
                         print.level = 2,
                         data = n.ml,
                         subset = n.ml$Treat == 1 & n.ml$indi_comp ==0
                         
)
summary(mnl.t1.other)

mnl.t2.indi_comp <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
               reflevel = 3,
               model = "mnl",  
               print.init = TRUE,
               print.level = 2,
               data = n.ml,
               subset = n.ml$Treat == 2 & n.ml$indi_comp ==1
               
)
summary(mnl.t2.indi_comp)

mnl.t2.other <- gmnl(formula = Choice ~ nobuy + price + gene  + color |0|0|0|0, 
                         reflevel = 3,
                         model = "mnl",  
                         print.init = TRUE,
                         print.level = 2,
                         data = n.ml,
                         subset = n.ml$Treat == 2 & n.ml$indi_comp ==0
                         
)
summary(mnl.t2.other)

####### table 2 #######
wtp.control.indi_comp <-coef(mnl.control.indi_comp)/(0-coef(mnl.control.indi_comp)[2])

start.control.indi_comp <- c(1, wtp.control.indi_comp[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.control.indi_comp <- gmnl(formula = Choice ~ price + nobuy + gene  + color|0|0|0|1,
                        reflevel = 3,
                        model="gmnl",
                        data=ml.wtp,
                        subset = ml.wtp$Treat==0 & ml.wtp$indi_comp ==1,
                        R=500,
                        ranp = c(gene = "n", color = "n"),
                        correlation = T,
                        fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                        panel = T, 
                        haltons = NA,
                        start = start.control.indi_comp,
                        print.init = 2,
                        print.level =2,
                        method = "bhhh",
                        rterlim = 500
                        
)
summary(wtp.rpl.control.indi_comp)
AIC(wtp.rpl.control.indi_comp)
BIC(wtp.rpl.control.indi_comp)
save(wtp.rpl.control.indi_comp, file = "wtp.rpl.control.indi_comp.Rdata")
load("wtp.rpl.control.indi_comp.Rdata")

se.cov.gmnl(wtp.rpl.control.indi_comp,sd=T)


wtp.t1.indi_comp <-coef(mnl.t1.indi_comp)/(0-coef(mnl.t1.indi_comp)[2])

start.t1.indi_comp <- c(1, wtp.t1.indi_comp[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t1.indi_comp <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                   reflevel = 3,
                   model="gmnl",
                   data=ml.wtp,
                   subset = ml.wtp$Treat==1 & ml.wtp$indi_comp ==1,
                   R=500,
                   ranp = c(gene = "n", color = "n"),
                   correlation = T,
                   fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                   panel = T, 
                   haltons = NA,
                   start = start.t1.indi_comp,
                   print.init = 2,
                   print.level =2,
                   method = "bhhh",
                   rterlim = 500
                   
)
summary(wtp.rpl.t1.indi_comp)
AIC(wtp.rpl.t1.indi_comp)
BIC(wtp.rpl.t1.indi_comp)
save(wtp.rpl.t1.indi_comp, file = "wtp.rpl.t1.indi_comp.Rdata")
load("wtp.rpl.t1.indi_comp.Rdata")

se.cov.gmnl(wtp.rpl.t1.indi_comp,sd=T)

wtp.t2.indi_comp <-coef(mnl.t2.indi_comp)/(0-coef(mnl.t2.indi_comp)[2])

start.t2.indi_comp <- c(1, wtp.t2.indi_comp[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t2.indi_comp <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                   reflevel = 3,
                   model="gmnl",
                   data=ml.wtp,
                   subset = ml.wtp$Treat==2 & ml.wtp$indi_comp ==1,
                   R=500,
                   ranp = c(gene = "n", color = "n"),
                   correlation = T,
                   fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                   panel = T, 
                   haltons = NA,
                   start = start.t2.indi_comp,
                   print.init = 2,
                   print.level =2,
                   method = "bhhh",
                   rterlim = 500
                   
)
summary(wtp.rpl.t2.indi_comp)
AIC(wtp.rpl.t2.indi_comp)
BIC(wtp.rpl.t2.indi_comp)
save(wtp.rpl.t2.indi_comp, file = "wtp.rpl.t2.indi_comp.Rdata")
load("wtp.rpl.t2.indi_comp.Rdata")

se.cov.gmnl(wtp.rpl.t2.indi_comp,sd=T)


wtp.control.other <-coef(mnl.control.other)/(0-coef(mnl.control.other)[2])

start.control.other <- c(1, wtp.control.other[c(1,3,4)],0.1,rep(0,3),0.1,0)
wtp.rpl.control.other <- gmnl(formula = Choice ~ price + nobuy + gene  + color|0|0|0|1,
                                  reflevel = 3,
                                  model="gmnl",
                                  data=ml.wtp,
                                  subset = ml.wtp$Treat==0 & ml.wtp$indi_comp ==0,
                                  R=500,
                                  ranp = c(gene = "n", color = "n"),
                                  correlation = T,
                                  fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                                  panel = T, 
                                  haltons = NA,
                                  start = start.control.other,
                                  print.init = 2,
                                  print.level =2,
                                  method = "bhhh",
                                  rterlim = 500
                                  
)
summary(wtp.rpl.control.other)
AIC(wtp.rpl.control.other)
BIC(wtp.rpl.control.other)
save(wtp.rpl.control.other, file = "wtp.rpl.control.other.Rdata")
load("wtp.rpl.control.other.Rdata")

se.cov.gmnl(wtp.rpl.control.other,sd=T)


wtp.t1.other <-coef(mnl.t1.other)/(0-coef(mnl.t1.other)[2])

start.t1.other <- c(1, wtp.t1.other[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t1.other <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                             reflevel = 3,
                             model="gmnl",
                             data=ml.wtp,
                             subset = ml.wtp$Treat==1 & ml.wtp$indi_comp ==0, 
                             R=500,
                             ranp = c(gene = "n", color = "n"),
                             correlation = T,
                             fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                             panel = T, 
                             haltons = NA,
                             start = start.t1.other,
                             print.init = 2,
                             print.level =2,
                             method = "bhhh",
                             rterlim = 500
                             
)
summary(wtp.rpl.t1.other)
AIC(wtp.rpl.t1.other)
BIC(wtp.rpl.t1.other)
save(wtp.rpl.t1.other, file = "wtp.rpl.t1.other.Rdata")
load("wtp.rpl.t1.other.Rdata")

se.cov.gmnl(wtp.rpl.t1.other,sd=T)

wtp.t2.other <-coef(mnl.t2.other)/(0-coef(mnl.t2.other)[2])

start.t2.other <- c(1, wtp.t2.other[c(1,3,4)],0.1,rep(0.1,3),0.1,0)
wtp.rpl.t2.other <- gmnl(formula = Choice ~ price + nobuy + gene  + color |0|0|0|1,
                             reflevel = 3,
                             model="gmnl",
                             data=ml.wtp,
                             subset = ml.wtp$Treat==2 & ml.wtp$indi_comp ==0,
                             R=500,
                             ranp = c(gene = "n", color = "n"),
                             correlation = T,
                             fixed= c(T, rep(F,3), F,rep(F,3), F,T),
                             panel = T, 
                             haltons = NA,
                             start = start.t2.other,
                             print.init = 2,
                             print.level =2,
                             method = "bhhh",
                             rterlim = 500
                             
)
summary(wtp.rpl.t2.other)
AIC(wtp.rpl.t2.other)
BIC(wtp.rpl.t2.other)
save(wtp.rpl.t2.other, file = "wtp.rpl.t2.other.Rdata")
load("wtp.rpl.t2.other.Rdata")

se.cov.gmnl(wtp.rpl.t2.other,sd=T)

##### table 3 #######

#between model
model1 <- c("wtp.rpl.control.indi_comp","wtp.rpl.control.indi_comp","wtp.rpl.t1.indi_comp")
model2 <- c("wtp.rpl.t1.indi_comp","wtp.rpl.t2.indi_comp","wtp.rpl.t2.indi_comp")

model1 <- c("wtp.rpl.control.other","wtp.rpl.control.other","wtp.rpl.t1.other")
model2 <- c("wtp.rpl.t1.other","wtp.rpl.t2.other","wtp.rpl.t2.other")

vari_compair <- c("gene", "color")
kr.poe(model1, model2, vari_compair)

# within model

model1 <- c("wtp.rpl.control.indi_comp","wtp.rpl.t1.indi_comp","wtp.rpl.t2.indi_comp")
model2 <- c("wtp.rpl.control.other","wtp.rpl.t1.other","wtp.rpl.t2.other")

vari_compair <- c("gene", "color")
kr.poe(model1, model2, vari_compair)


